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Preface 



In my lectures at the Les Houches Summer School 2008, I discussed central concepts 
of computational statistical physics, which I felt would be accessible to the very cross- 
cultural audience at the school. 

I started with a discussion of sampling, which lies at the heart of the Monte Carlo 
approach. I specially emphasized the concept of perfect sampling, which offers a syn- 
thesis of the traditional direct and Markov-chain sampling approaches. The second 
lecture concerned classical hard-sphere systems, which illuminate the foundations of 
statistical mechanics, but also illustrate the curious difficulties that beset even the 
most recent simulations. I then moved on, in the third lecture, to quantum Monte 
Carlo methods, that undcrly much of the modern work in bosonic systems. Quantum 
Monte Carlo is an intricate subject. Yet one can discuss it in simplified settings (the 
single-particle free propagator, ideal bosons) and write direct-sampling algorithms for 
the two cases in two or three dozen lines of code only. These idealized algorithms illus- 
trate many of the crucial ideas in the field. The fourth lecture attempted to illustrate 
aspects of the unity of physics as realized in the Ising model simulations of recent 
years. 

More details on what I discussed in Les Houches, and wrote up (and somewhat 
rearranged) here, can be found in my book, "Statistical Mechanics: Algorithms and 
Computations" (SMAC), as well as in recent papers. Computer programs are available 
for download and perusal at the book's web site www.smac.lps.ens.fr. 
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1 

Sampling 



1.1 Direct sampling, sample transformation 

As an illustration of what is meant by sampling, and how it relates to integration, we 
consider the Gaussian integral: 

' exp (-x 2 /2) = 1. (1.1) 
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This integral can be computed by taking its square: 
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and then switching to polar coordinates (dxdy = rdrd</>), 
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as well as performing the substitutions r /2 = T (rdr = dT) and exp (— T) = $ 
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In our context, it is less important that we can do the integrals in eqn p.6[) analytically, 
than that <f> and \P can be sampled as uniform random variables in the interval [0, 27r] 
(for </>) and in [0, 1] (for <3>). Samples = ran(0, 2-k) and 'J = ran(0, 1), are readily 
obtained from the random number generator ran (a, b) lingering on any computer. We 
can plug these random numbers into the substitution formulas which took us from 
eqn (|1.2|) to eqn (|1.6p . and that take us now from two uniform random numbers <f> 
and VP to Gaussian random numbers x and y. We may thus apply the integral trans- 
formations in the above equation to the samples, in other words perform a "sample 
transformation" . This is a practical procedure for generating Gaussian random num- 
bers from uniform random numbers, and we best discuss it as what it is, namely an 
algorithm. 
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procedure gauss 

ran(0,27r) 

ran (0,1) 

-log \t 
/2T 

cos 
y <— rsin 
output {x, y} 

Algorithm 1 gauss. SMAC pseudocode program for transforming two uniform ran- 
dom numbers into two independent uniform Gaussian random numbers x,y. 
Except for typography, SMAC pseudocode is close to the Python programming lan- 
guage. 

SMAC pseudocode can be implemented in many computer languages 0. Of partic- 
ular interest is the computer language Python, which resembles pseudocode, but is 
executable on a computer exactly as written. We continue in these lectures using and 
showing Python code, as in algorithm 1.1. 

Algorithm 1.1 gausstest.py 



from random import uniform as ran , gauss 
from math import sin , cos , sqrt , log , pi 
def gausstest ( sigma ) : 

phi = ran(0,2*pi) 

Upsilon = — log(ran(0. ,1-)) 

r = sigma* sqrt (2* Upsilon ) 

x = r*cos(phi) 

y = r * sin ( phi ) 

return x,y 
print gausstest ( 1 . ) 



Direct-sampling algorithms exist for arbitrary one-dimensional distributions (see 
SMAC Sect. 1.2.4). Furthermore, arbitrary discrete distributions {7ri, . . . ,ttk} can be 
directly sampled, after an initial effort of about K operations with ~ log 2 K operation 
by "Tower sampling" (see SMAC Sect. 1.2.3). This means that sampling a distribu- 
tion made up of, say, one billion terms takes only about 30 steps. We will use this 
fact in Section 13 . 31 for the direct sampling algorithm for ideal bosons. Many trivial 
multi-dimensional distribution (as for example non-interacting particles) can also be 
sampled. 

Direct-sampling algorithms also solve much less trivial problems as for example 
the free path integral, ideal bosons, and the two-dimensional Ising model, in fact, 
many problems which possess an exact analytic solution. These direct-sampling algo- 
rithms are often the race-car engines inside general-purpose Markov-chain algorithms 
for complicated interacting problems. 

Let us discuss the computation of integrals derived from the sampled ones Z = 
J dx 7r(x) (in our example, the distribution ir(x) is the Gaussian) 

lr The SMAC web page www.smac.lps.ens.fr provides programs in language ranging from Python, 
Fortran, C, and Mathematica to TI basic, the language of some pocket calculators. 
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JdxO(x)n(x) ]_^ n( , 
fdx7r(x) ~N^ U{Xlh 



where the points Xi are sampled from the distribution tt(x). This approach, as shown, 
allows to compute mean values of observables for a given distribution n. One can also 
bias the distribution function it using Markov-chain approaches. This gives access to 
a large class of very non-trivial distributions. 

1.2 Markov-chain sampling 

Before taking up the discussion of elaborate systems (hard spheres, bosons, spin 
glasses), we first concentrate on a single particle in a finite one-dimensional lat- 
tice k = 1, . . . ,N (sec Fig. II. ip . This case is even simpler than the aforementioned 
Gaussian, because the space is discrete rather than continuous, and because the site- 
occupation probabilities {tti, . . . , ttn} are all the same 

TT k = — Vfc. 

We can sample this trivial distribution by picking k as a random integer between 1 
and N. Let us nevertheless study a Markov-chain algorithm, whose diffusive dynam- 
ics converges towards the probability distribution of eqn (|1.2p . For concreteness, wc 
consider the algorithm where at all integer times tg]- oo, oof, the particle hops with 
probability | from one site to each of its neighbors: 

p k ->k+i =Pk^k-i = 1/3 (if possible). (1.7) 

The probabilities to remain on the site are 1/3 in the interior and, at the boundaries, 
Pi^i = Pn^n = 2/3. 
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Fig. 1.1 Markov-chain algorithm on a five-site lattice. A single particle hops towards 
a site's neighbor with probability 1/3, and remains on the site with probability 1/3 
(probability 2/3 on the boundary). 

These transition probabilities satisfy the notorious detailed balance condition 

KkPk^i = npi->k (1-8) 

for the constant probability distribution ixk. Together with the crgodicity of the al- 
gorithm this guarantees that in the infinite-time limit, the probability itk to find the 
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particle at site k is indeed independent of k. With appropriate transition probabilities, 
the Metropolis algorithm would allow us to sample any generic distribution 7Tfc or 7r(x) 
in one or more dimensions (sec SMAC Sect. 1.1.5). 

Markov-chain methods arc more general than direct approaches, but the price to 
pay is that they converge to the target distribution tt only in the infinite-time limit. 
The nature of this convergence can be analyzed by the transfer matrix T 1,1 (the matrix 
of transition probabilities, in our case a 5 x 5 matrix): 



T 



{Pk^i} 



(2 1 0\ 
1110 
1110 
111 

\0 1 2) 



(1.9) 



The eigenvalues of the transfer matrix T are (1, \ ± V5/6, \ ± \/5/6). The largest 
eigenvalue, equal to one, expresses the conservation of the sum of all probabilities 
^(■7ri + • • • + ttq) = 1. Its corresponding eigenvector is | • o o oo) + • • • + | o o o o»), by 
construction, because of the detailed-balance condition eqn (|4.2[) . The second-largest 
eigenvalue, A2' 1 = \ + V5/6 = 0.8727 governs the decay of correlation functions 
at large times. This is easily seen by computing the probability vector 7r(to + t) = 
{7ri(t), 775(4)} which can be written in terms of the eigenvalues and eigenvectors 

Let us suppose that the simulations always star10 on site k = 1, and let us decom- 
pose this initial configuration onto the eigenvectors tt(0) = a\Tx\ + ■ • ■ + aNW e N . 

n(t) = (T 1 * 1 )* tt(0) = + a % \\^\ + ... 

At large times, corrections to the equilibrium state vanish as A£, so that site proba- 
bilities approach the equilibrium value as exp (r/r corr ) with a time constant 



W = l/|logA a | ^7.874 



(1.10) 



(see SMAC Sect. 1.1.4, p. 19f). 

We retain that the exponential convergence of a Monte Carlo algorithm is charac- 
terized by a scale, the convergence time r corr . We also retain that convergence takes 
place after a few r cor r, in our example say 3 x 7.87 ~ 25 iterations (there is absolutely 
no need to simulate for an infinite time). If our Markov chains always start at t = on 
the site 1, it is clear that for all times t' < 00, the occupation probability TTi(t') of site 
1 is larger than, say, the probability TTs(t') to be on site 5. This would make us believe 
that Markov chain simulations never completely decorrelate from the initial condition, 
and are always somehow less good than direct-sampling. This belief is wrong, as we 
shall discuss in the next section. 



1.3 Perfect sampling 

We need to simulate for no more than a few T cor r, but we must not stop our calculation 
short of this time. This is the critical problem in many real-life situations, where we 

2 we should assume that the initial configuration is different from the stationary solution because 
otherwise all the coefficients <X2,...,atf would be zero. 



Perfect sampling 5 



cannot compute the correlation time as reliably as in our five-site problem: r corr may 
be much larger than we suspect, because the empirical approaches for determining 
correlation times may have failed (sec SMAC Sect, f .3.5). In contrast, the problem of 
the exponentially small tail for r S> t co1t is totally irrelevant. 

Let us take up again our five-site problem, with the goal of obtaining rigorous 
information about convergence from within the simulation. As illustrated in Fig. 11.2) 
the Markov-chain algorithm can be formulated in terms of time sequences of random 
maps: Instead of prescribing one move per time step, as we did in Fig. 11.11 we now 
sample moves independently for all sites k and each time t. At time to, for example, 
the particle should move straight from sites 1, 2 and 5 and down from sites 3 and 4, 
etc. Evidently, for a single particle, there is no difference between the two formulations, 
and detailed balance is satisfied either way. With the formulation in terms of random 
maps, we can verify that from time to + ^coup on, all initial conditions generate the 
same output. This so-called "coupling" is of great interest because after the coupling 
time T C oup, the influence of the initial condition has completely disappeared. In the rest 
of this lecture, we consider extended Monte Carlo simulations as the one in Fig. 11.21 
with arrows drawn for each site and time. 



f f +x coup 

Fig. 1.2 Extended Monte Carlo simulation on N = 5 sites. In this example, the 
coupling time is r coup = 11. 



The coupling time r coup is a random variable (r coup = 11 in Fig. II. 2[) whose dis- 
tribution 7r(r COU p) vanishes exponentially in the limit r coup — » oo because the random 
maps at different times are independent. 

The extended Monte Carlo dynamics describes a physical system with from 1 to 
N particles, and transition probabilities, from eqn (|1.7p . as for example: 



p fw (| o o • o.) -> I o • o «o)) = 1/9 
p fw (| oo.o.) -> I oo.o.)) = 2/9 

, 71 " ' (1.11) 

p fw (| oo.o.) -> 000.0)) = 1/9 
p fw (| o o •••)-> I o o o .o)) = 1/27, 



etc., We may analyze the convergence of this system through its transfer matrix T 
(in our case, a 31 x 31 matrix, because of the 31 non-empty states on five sites): 
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Apl,! y2,l 

T 2 > 2 T 3 > 2 
T 3 ' 3 

\ 



where the block T k ' 1 connects states with k particles at time t to states with / < k 
particles at time t + 1. The matrix T fw has an eigenvector with eigenvalue 1 which 
describes the equilibrium, and a second-largest eigenvalue which describes the ap- 
proach to equilibrium, and yields the coupling time r coup . Because of the block- 
triangular structure of T fw , the eigenvalues of T x > x constitute the single-particle sector 
of T fw , including the largest eigenvalue A( w = 1, with corresponding left eigenvector 
| • o o oo) + ■ • • + 1 o o o o») . The second-largest eigenvalue of T belongs to the (2, 2) sec- 
tor. It is given by AJ, W = 0.89720 > A% , and it describes the behavior of the coupling 
probability 7t(t coup ) for large times. 

Markov chains shake off all memory of their initial conditions at time r coup , but 
this time changes from simulation to simulation: it is a random variable. Ensemble 
averages over (extended) simulations starting at to, and ending at <o + r thus contain 
mixed averages over coupled and "not yet coupled" runs, and only the latter carry 
correlations with the initial condition. To reach pure ensemble averages over coupled 
simulations only, one has to start the simulations at time t = — oo, and go up to t = 0. 
This procedure, termed "coupling from the past" (Propp and Wilson 1996), is familiar 
to theoretical physicists because in dynamical calculations, the initial condition is often 
put to — oo. Let us now see how this trick can be put to use in practical calculations. 
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Fig. 1.3 Extended simulation on iV = 5 sites. The outcome of this infinite simulation, 
from t — —oo through t = 0, is k — 2. 



An extended simulation "from the past" is shown in Fig. 11.31 It has run for an 
infinite time so that all the eigenvalues of the transfer matrix but Ai = 1 have died 
away and only the equilibrium state has survived. The configuration at t = is thus 
a perfect sample, and each value k (fee {1, . . . , 5} is equally likely, if we average over 
all extended Monte Carlo simulations. For the specific simulation (choice of arrows in 
Fig. |1.3|) . we know that the simulation must pass by one of the five points at time 
t = to = —10). However, the Markov chain couples between t — to and t = 0, so that 
we know that k(t = 0) = 2. If it did not couple at t = to, we would have to provide 
more information (draw more arrows, for t = to — 1, to — 2, . . . ) until the chain couples. 
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Up to now, we have considered the forward transfer matrix, but there is also a 
backward matrix T bw , which describes the propagation of configurations from t = 
back to t = -1, t = -2, etc. A hw is similar to the matrix A lw ~ A hw ( the two 
matrices have identical eigenvalues), because the probability distribution to couple 
between time t = — \to\ and t = equals the probability to couple between time t = 
and t = +|to|. This implies that the distributions of coupling times in the forward and 
backward directions are identical. 

We have also seen that the eigenvalue corresponding to the coupling time is larger 
than the one for the correlation time. For our one-dimensional diffusion problem, this 
can be proven exactly. More generally, this is because connected correlation functions 
decay as 

(O(0)O(t)) e ~ e"^— . 

Only the non-coupled (extended) simulations contribute to this correlation function, 
and their proportion is equal to cxp (— t/r coup ). We arrive at 



exp [-* (l/r COU p - l/r corr )] ~ (O(O)O(t)) 



non-coup 



and T coup > T corr because even the non-coupling correlation functions should decay 
with time. 

Finally, we have computed in Fig. 11.21 and in Fig. 11.31 the coupling time by follow- 
ing all possible initial conditions. This approach is unpractical for more complicated 
problems, such as the Ising model on N sites with its 2 N configurations. Let us show in 
the five-site model how this problem can sometimes be overcome: It suffices to define 
a new Monte Carlo dynamics in which trajectories cannot cross, by updating, say, 
even lattice sites (k = 2,4) at every other time to, to + 2, . . . and odd lattice sites 
(k = 1, 3, 5) at times to + 1, to + 3, . . . . The coupling of the trajectories starting at sites 
1 and 5 obviously determines r coup . 
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Fig. 1.4 Extended Monte Carlo simulation with odd and even times. Trajectories 
cannot cross, and the coupling of the two extremal configurations (with A: (to) = 1 and 
k(to = 5)) determines the coupling time. 

In the Ising model, the ordering relation of Fig. [L4] survives in the form of a "half- 
order" between configurations (see SMAC Sect. 5.2.2), but in the truly interesting 
models, such as spin glasses, no such trick is likely to exist. One must really supervise 
the 2 N configurations at t = to. This non-trivial task has been studied extensively (see 
Chanal and Krauth 2008). 
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2.1 Molecular dynamics 

Before statistical mechanics, not so long ago, there only was classical mechanics. The 
junction between the two has fascinated generations of mathematicians and physicists, 
and nowhere can it be studied better than in the hard-sphere system: N particles in a 
box, moving about like billiard balls, with no other interaction than the hard-sphere 
exclusion (without friction or angular momentum) . For more than a century, the hard- 
sphere model has been a prime example for statistical mechanics and a parade ground 
for rigorous mathematical physics. For more than fifty years, it has served as a test 
bed for computational physics, and it continues to play this role. 

For concreteness, we first consider four disks (two-dimensional spheres) in a square 
box with walls (no periodic boundary conditions), as in Fig. 12.11 From an initial 
configuration, as the one at t = 0, particles fly apart on straight trajectories either 
until one of them hits a wall or until two disks undergo an elastic collision. The 
time for this next "event" can be computed exactly by going over all disks (taken by 
themselves, in the box, to compute the time for the next wall collisions) and all pairs 
of disks (isolated from the rest of the system, and from the walls, to determine the 
next pair collisions), and then treating the event closest in time (see SMAC Sect. 2.1). 

wall collision pair collision 
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Fig. 2.1 Event-driven molecular dynamics simulation with four hard disks in a square 
box. 

The event-chain algorithm can be implemented in a few dozen lines, just a few too 
many for a free afternoon in Les Houches (program listings are available on the SMAC 
web site). It implements the entire dynamics of the iV-particle system without time 
discretization, because there is no differential equation to be solved. The only error 
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committed stems from the finite-precision arithmetic implemented on a computer. 

This error manifests itself surprisingly quickly, even in our simulation of four disks 
in a square: typically, calculations done in 64bit precision (15 significant digits) get out 
of step with other calculations from identical initial conditions with 32bit calculations 
after a few dozen pair collisions. This is the manifestation of chaos in the hard-sphere 
system. The appearance of numerical uncertainties for given initial conditions can be 
delayed by using even higher precision arithmetic, but it is out of the question to 
control a calculation that has run for a few minutes on our laptop, and gone through 
a few billion collisions. Chaos in the hard sphere model has its origin in the negative 
curvature of the sphere surfaces, which magnifies tiny differences in the trajectory at 
each pair collision and causes serious roundoff errors in computation^. 

The mathematically rigorous analysis of chaos in the hard-disks system has met 
with resounding success: Sinai (1970) (for two disks) and Simanyi (2003, 2004) (for 
general hard-disk and hard-sphere systems) were able to mathematically prove the 
foundations of statistical physics for the system at hand, namely the equiprobability 
principc for hard spheres: This means that during an infinitely long molecular dynamics 
simulation, the probability density satisfies the following: 

{probability of configuration with 1 , . , . 

r ,i I r , - 1 f OC7r(Xi,...,Xjv)dXi,...,dXjv, 

[xi,xi + dxij, . . . , [xat,xjv + dxjvj J 

where 

J 1 if configuration legal (|xfc — xj| > 2cr for k ^ I) 
7r(xi,...,xjv) = < . (2.1) 

I otherwise 

An analogous property has been proven for the velocities: 

fl if velocities legal (£)v| = E kin / (2m)) 
tt(vi, . . .,v N ) = < . {2.2) 

I otherwise 

Velocities are legal if they add up to the correct value of the conserved kinetic energy, 
and their distribution is constant on the surface of the 2N dimensional hypersphcrc 
of radius y/E^J(2m) (for disks). 

The two above equations contain all of equilibrium statistical physics in a nutshell. 
The equal-probability principle of eqn (|2.ip relates to the principle that two config- 
urations of the same energy have the same statistical weight. The sampling problem 
for velocities, in eqn (|2.2j) can be solved with Gaussians, as discussed in SMAC Sect. 
1.2.6. It reduces to the Maxwell distribution for large N (sec SMAC Sect. 2.2.4) and 
implies the Boltzmann distribution (see SMAC Sect. 2.3.2). 

2.2 Direct sampling, Markov chain sampling 

We now move on from molecular dynamics simulations to the Monte Carlo sampling. 
To sample N disks with the constant probability distribution of eqn (|2.1[) . we uniformly 

1 as it causes humiliating experiences at the billiard table 
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Algorithm 2.1 direct-disks. py 



1 


from random import uniform as ran 




2 


sigma = 0.20 




3 


condition = False 




4 


while condition = False: 




5 


L = [( ran ( sigma , 1 — sigma ), ran ( sig 


ma , 1 — sigma ) ) ] 


6 


for k in range ( 1 , 4) : #4 P artic 


les considered 


7 


b = ( ran ( sigma , 1 — sigma ), ran ( sig 


ma , 1 — sigma ) ) 


8 


min_dist = min ( ( b [0] — x [0] ) * * 2 + 


(b[l]-x[l])**2 for x in L) 


9 


if min_dist < 4*sigma**2: 




10 


condition = False 




11 


break 




12 


else : 




13 


L . append ( b ) 




14 


condition = True 




IS 


print L 





throw a set of N particle positions {xi, . . . , xjy} into the square. Each of these sets of 
N positions is generated with equal probability. We then sort out all those sets that 
are no legal hard-sphere configurations. The remaining ones (the gray configurations 
in Fig. 12. 2|) still have equal probabilities, exactly as called for in eqn (|4.ip . In the 
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Fig. 2.2 Direct-sampling Monte Carlo algorithm for hard disks in a box without 
periodic boundary conditions (see Alg. direct-disks .py). 



Python programming language, we can implement this algorithm in a few lines (see 
Alg. direct-disks .py): one places up to N particles at random positions (see line 7 of 
Alg. direct-disks .py). If two disks overlap, one breaks out of this construction and 
restarts with an empty configuration. The rejection rate p aC cept{N, a) of this algorithm 
(the probability to generate legal (gray) configurations in Fig. 12. 2p : 

number of valid configurations with radius a ^ 
accept number of valid configurations with radius a 

is exponentially small both in the particle number N and in the density of particles, 
and this for physical reasons (see the discussion in SMAC Sect. 2.2.2)). 

For the hard-sphere system, Markov-chain methods are much more widely ap- 
plicable than the direct-sampling algorithm. In order to satisfy the detailed balance 
condition of cqn (|4.2j) , we must impose that the probability to move from a legal con- 
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Algorithm 2.2 markov-disks.py 

1 from random import uniform as ran , choice 

2 L = [(0. 25, 0.25), (0.75, 0.25), (0.25, 0.75), (0.75, 0.75)] 

3 sigma , delta =0.20 ,0. 15 

4 for iter in range ( 1000000) : 

5 a = choice (L) 

6 L . remove (a) 

7 b = ( a [ ] + ran( — delta, delta), a [ 1 ] + ran(— delta , delta)) 

8 min.dist = min ( ( b [0] -x [0] ) * * 2 + ( b [1] - x [ 1 ] ) * * 2 for x in L) 

9 box.cond = min(b[0] ,b[l]) < sigma or max(b[0] ,b[l]) >1— sigma 

10 if box.cond or min_dist < 4*sigma**2: 

11 L . append (a) 

12 else : 

13 L . append (b) 

14 print L 



figuration a to another, b, must be the same as the probability to move from b back to a 
(see Fig. l2.3|) . This is realized most easily by picking a random disk and moving it inside 
a small square around its original position, as implemented in Alg. markov-disks .py. 
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Fig. 2.3 Markov-chain Monte Carlo algorithm for hard disks in a box without periodic 
boundary conditions. 



The Monte-Carlo dynamics of the Markov-chain algorithm, in Fig. 12.31 superficially 
resembles the molecular-dynamics in Fig. 12.11 Let us work out the essential differences 
between the two: In fact, the Monte Carlo dynamics is diffusive: It can be described 
in terms of diffusion constants and transfer matrices, and convergence towards the 
equilibrium distribution of eqn (|2.1| is exponential. In contrast, the dynamics of the 
event-chain algorithm is hydrodynamic (it contains eddies, turbulence, etc, and their 
characteristic timescales). Although it converges to the same equilibrium state, as dis- 
cussed before, this convergence is algebraic. This has dramatic consequences, especially 
in two dimensions (sec SMAC Sect. 2.2.5), even though the long-time dynamics in a 
finite box is more or less equivalent in the two cases. This is the fascinating subject 
of long-time tails, discovered by Alder and Wainwright (1970), which for lack of time 
could not be covered in my lectures (sec SMAC Sect. 2.2.5). 

The local Markov-chain Monte Carlo algorithm runs into serious trouble at high 
density, where the Markov chain of configurations effectively gets stuck during long 
times (although it remains ergodic). The cleanest illustration of this fact is obtained 
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disk k same disk 




z = ... i = 25600000000 



Fig. 2.4 Initial and final configurations for 256 disks at density r\ = 0.72 (much 
smaller than the close-packing density) after 25.6 billion iterations. The sample at 
iteration i = 25.6 x 10 9 remembers the initial configuration. 

by starting the run with an easily recognizable initial configuration, as the one shown 
in Fig. 12.41 which is slightly tilted with respect to the x-axis. We see in this example 
that T corr 3> 2.56 x 10 9 iterations, even though it remains finite and can be measured 
in a Monte Carlo calculation. 

Presently, no essentially faster algorithm than the local algorithm is known for 
uniform hard spheres, and many practical calculations (done with very large numbers 
of particles at high density) are clearly un-converged. Let us notice that the slowdown 
of the Monte Carlo calculation at high density has a well-defined physical origin: the 
slowdown of the single-particle diffusion at high density. However, this does not exclude 
the possibility of much faster algorithms, as we will discuss in the fourth lecture. 

2.3 Cluster algorithms, birth-and-death processes 

In the present section, we explore Monte Carlo algorithms that are not inspired by the 
physical process of single-particle diffusion underlying the local Markov-chain Monte 
Carlo algorithm. Instead of moving one particle after the other, cluster algorithms 
construct coordinated moves of several particles at a time, from one configuration, a, to 
a very different configuration, b in one deterministic step. The pivot cluster algorithm 
(Dress and Krauth 1995, Krauth and Moessner 2003) is the simplest representative of 
a whole class of algorithms. 




a a' a" ... b 

Fig. 2.5 A cluster move about a pivot axis involving five disks. The pivot-axis posi- 
tion and its orientation (horizontal, vertical, diagonal) are chosen randomly. Periodic 
boundary conditions allow us to scroll the pivot into the center of the box. 
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Algorithm 2.3 pocket-disks. py in a periodic box of size lxl 



from random import uniform as ran , choice 

import box_it , dist # see SMAC web site for periodic distance 
Others = [(0.2 5 ,0.25), (0.25, 0.75), (0.75, 0.25), (0.75, 0.75)] 
sigma.sq =0.15**2 
for iter in range ( 1 00 00 ) : 

a = choice ( Others ) 

Others . remove ( a ) 

Pocket = [a] 

Pivot = (ran (0 , 1) , ran (0,1)) 
while Pocket != [] : 

a = Pocket, pop () 

a = T(a , Pivot ) 

for b in Others [:] : # "Others [:]" is a copy of "Others" 
if dist(a,b) < 4*sigma_sq: 

Others . remove (b) 

Pocket . append (b) 
Others . append (a) 
print Others 



In this algorithm, one uses a "pocket" , a stack of disks that eventually have to 
move. Initially, the pocket contains a random disk. As long as there are disks in the 
pocket, one takes one of them out of it and moves it. It gets permanently placed, and all 
particles it overlaps with are added to the pocket (see Fig. 12.51 the "pocket particles" 
are colored in dark). In order to satisfy detailed balance, the move a — > b = T(a) 
must have a Z2 symmetry a = T 2 (a) (as a reflection around a point, an axis, or a 
hypcrplane), such that moving it twice brings each particle back to its original position 
(see SMAC Sect. 2.5.2). In addition, the transformation T must map the simulation box 
onto itself. For example, we can use a reflection around a diagonal in a square or cube 
box, but not in a rectangle or cuboid. Periodic boundary conditions are an essential 
ingredient in this algorithm. In Python, the pocket algorithm can be implemented in 
a dozen lines of code (see Alg. pocket-disks .py). 




b 



Fig. 2.6 Single iteration (from a to b) of the cluster algorithm for a binary mixture of 
hard spheres. The symmetry operation is a flip around the y axis shown. Transformed 
disks are marked with a dot (see Buhot and Krauth, 1999). 



The pivot-cluster algorithm fails at high density, say, at the condition of Fig. 12. A\ 
where the transformation T simply transforms all particles in the system. In that case 



14 Classical hard-sphere systems 



entire system without changing the relative particle positions. However, the algorithm 
is extremely efficient for simulations of monomcr-dimer models or in binary mixtures, 
among others. An example of this is given in Fig. 12.61 

At the end of this lecture, let us formulate hard-sphere systems with a grand- 
canonical partition function and fugacity A: 

oo 

Z = ^2 \ N ir(xi, . . .,x N ). 
n=o 

and discuss the related Markov-chain Monte Carlo algorithms in terms of birth-and- 
death processes (no existential connotation intended): Between two connected con- 
figurations, the configuration can only change through the appearance ( "birth" ) or 
disappearance ("death") of a disk (see Fig. 12. 7j) It follows from the detailed balance 



oo 




oo 




o 


oo 




oo 




oo 



a move b 



Fig. 2.7 Birth-and-death process for grand-canonical hard disks. 

condition eqn (|4.2|) that the probability to ir(a) = A-7r(6), so that p(b — * a) = \p(a — > b). 
This means that to sample eqn (|2.3|) . we simply have to install one "death" probability 
(per time interval dt) for any particle, as well as a birth probability for creating a new 
particle anywhere in the system (this particle is rejected if it generates an overlap). 

As in the so-called "fastcr-than-the-clock" algorithms (sec SMAC Sect. 7.1) one 
does not discretize time r — > A T but rather samples lifetimes and birth times from 
their exponential distributions. In Fig. 12.81 we show a time-space diagram of all the 
events that can happen in one extended simulation (as in the first lecture), in the 
coupling-from-the-past framework, as used by Kendall and Moller (2000). We do not 
know the configuration of disks but, on a closer look, we can deduce it, starting from 
the time to indicated. 




(now) 

Fig. 2.8 One-dimensional time-space diagram of birth-and-death processes for hard 
disks in one dimension. 



3 

Quantum Monte Carlo simulations 



3.1 Density matrices, naive quantum simulations 

After the connection between classical mechanics and (classical) statistical mechanics, 
we now investigate the junction between statistical mechanics and quantum physics. 
We first consider a single particle of mass m in a harmonic potential V(x) = \(jJX 2 , 
for which wave functions and energy eigenvalues arc all known (see Fig. 13.11 we use 
units H = m = u> = I). 




1 i — 1 r 

-5 5 



position x 

Fig. 3.1 Harmonic-oscillator wave functions ip n , shifted by the energy eigenvalue E n . 

As before, we are interested in the probability tt(x) for a particle to be at position 
x (compare with eqn (|2.ip ). This probability can be assembled from the statistical 
weight of level k, ir(Ek) oc cxp (—(3Ek) and from the quantum mechanical probability 
tjjk(x)ipk(x) to be at position x while in level k, 

<x) = -S- Y, ex P irPEk) *k(x)n(x) ■ (3-1) 



This probability involves a diagonal clement of the density matrix given by p(x, x' , (3) = 
J2k ex P ^ k(x)^* k (x) , whose trace is the partition function: 
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Z(p)= I dx p(x,x,0)=J2 ex P{-P E k) 



(see SMAC Sect 3.1.1). For the free particle and the harmonic oscillator, we can indeed 
compute the density matrix exactly (see the later eqn (|3.5p ) but in general, eigenstates 
or energies are out of reach for more complicated Hamiltonians. The path integral ap- 
proach obtains the density matrix without knowing the spectrum of the Hamiltonian, 
by using the convolution property 

p any {x, x',f3) = [ dx" p any (x, x", (3') p any (x", x', (3 - fi) (convolution), (3.2) 



which yields the density matrix at a given temperature through a product of density 
matrices at higher temperatures. In the high-temperature limit, the density matrix for 
the Hamiltonian H = H bee + V is given by the Trotter formula 

p any (x,x',(3) y-^e-i^^p^ix,^^) e -hPV{x') (Trotter formula). (3.3) 

(compare with SMAC Sect. 3.1.2). For concreteness, we continue our discussion with 
the harmonic potential V = \x 2 , although the discussed methods are completely 
general. Using eqn (|3.2[) repeatedly, one arrives at the path-integral representation of 
the partition function in terms of high-temperature density matrices: 



Z{0) — J dxp(x,x,0)= j dxo ••• / dxN-i 

p(x ,xi,A T ) . . .p(x N -i,x ,A T ) (3.4) 
y * 

7t(K0,...,XW-i) 

where A T = f3/N . In the remainder of this lecture we will focus on this multiple in- 
tegral but we are again more interested in sampling (that is, in generating "paths" 
{xo, . . . ,Xjv_i} with probability ir(xo, . . . ,xn-i)) than in actually computing it. A 
naive Monte Carlo sampling algorithm is set up in a few lines of code (see the 
Alg. naive-harmonic-path.py, on the SMAC web site). In the integrand of eqn (|3.4p , 
one chooses a random point k £ [0, N — 1] and a uniform random displacement 
5 X e (compare with Fig. 13. 2[) . The acceptance probability of the move de- 

pends on the weights ir(xo, ■ ■ ■ , xjv-i), thus both on the free density matrix part, and 
on the interaction potential. 

The naive algorithm is extremely slow because it moves only a single "bead" k 
out of N, and not very far from its neighbors k — 1 and k + 1. At the same time, 
displacements of several beads at the same time would rarely be accepted. In order to 
go faster through configuration space, our proposed moves must learn about quantum 
mechanics. This is what we will teach them in the following section. 

3.2 Direct sampling of a quantum particle: Levy construction 

The problem with naive path-sampling is that the algorithm lacks insight: the prob- 
ability distribution of the proposed moves contains no information about quantum 
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; = 



a a (+ move) b 

Fig. 3.2 Move of the path {xq, • ■ • , xn-i} via element k G [0, N — 1], sampling the 
partition function in eqn (|3.4p . For k — 0, both and N are displaced in the same 
direction. The external potential is indicated. 

mechanics. However, this problem can be solved completely for a free quantum parti- 
cle and also for a particle in a harmonic potential V(x) = ^x 2 , because in both cases 
the exact density matrix for a particle of mass m = 1, with h= 1, is a Gaussian: 



p(x,x',j3) 



i 



exp 



{x-x'f 



2/3 



(free particle) 



^/2n sinli 3 



exp 



(x+x' 



tanh § — 



_ (x-x') 2 



cothf 



(osc.) 



(3.5) 



The distribution of an intermediate point, for example in the left panel of Fig. 1373 
is given by 



7r(xi|x ,a;Ar) 



P\ Xi,X N , 



Gaussian, see eqn ]3.5| l Gaussian, 



(N-1)0 
N 

see eqn I l3.5t 



As a product of two Gaussians, this is again a Gaussian, and it can be sampled directly. 
After sampling x\, one can go on to X2, etc., until the the whole path is constructed, 
and without rejections (Levy 1940). 




k= 1 




k = 2 



k=3 




k = 4 



1 




k = 5 




output 



Fig. 3.3 Levy construction of a path contributing to the free-particle density matrix 
p(xo, xe, /3). The intermediate distributions, all Gaussians, can be sampled directly. 
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Algorithm 3.1 levy-free-path. py 




1 


from math import pi , exp , sqrt 




2 


from random import gauss 




3 


N = 10000 




4 


beta = 1. 




5 


Dcl.tau = bcta/N 




6 


x = [0 for k in range (N+l)] 




7 


y = range (N+l) 




8 


for k in range (1,N): 




9 


Dcl.p = (N-k)* Del.tau 




10 


x_mean = ( Del_p *x [ k — 1] + Del.tau *x [N] ) / ( Del.tau - 


- Del.p) 


11 


sigma = sqrt ( 1 ./( 1 ./ Del.tau + l./Del_p)) 




112 


x[k] = gauss ( x.mean , sigma ) 




13 


print x 





The Levy construction is exact for a free particle in a harmonic potential and of 
course also for free particles (as shown in the Python code). Direct-sampling algorithms 
can be tailored to specific situations, such as periodic boundary conditions or hard 
walls (see SMAC Sections 3.3.2 and 3.3.3). Moreover, even for generic interactions, the 
Levy construction allows to pull out, treat exactly, and sample without rejections the 
free-particle Hamiltonian. In the generic Hamiltonian H = Hq + V , the Metropolis 
rejection then only takes care of the interaction term V. Much larger moves than before 
become possible and the phase space is run through more efficiently than in the naive 
algorithm, where the Metropolis rejection concerned the entire H. This makes possible 
nontrivial simulations with a very large number of interacting particles (see Krauth 
1996, Holzmann and Krauth 2008). We will illustrate this point in the following section, 
showing that iV-body simulations of ideal (non-interacting) bosons (of arbitrary size) 
can be done by direct sampling without any rejection. 

We note that the density matrices in eqn (|3.5| . for a single particle in a d-dimcnsional 
harmonic oscillator yield the partition functions 

z h -°-((3) = Jdxp(x,x,p) = [1 -exp(-^)r d (with£o=0) (3.6) 

where we use the lowercase symbol z((3) in order to differentiate the one-particle 
partition function from its A^-body counterpart Z(Ji) which we will need in the next 
section. 

3.3 Ideal bosons: Landsberg recursion and direct sampling 

The density matrix for N distinguishable particles p dlst ({xi, . . . , xn}, {x[, . . . , x' N }, 0) 
is assembled from normalized TV-particle wavefunctions and their associated energies, 
as in eqn (|3.ip . The density matrix for indistinguishable particles is then obtained 
by symmetrizing the density matrix for N distinguishable particles. This can be done 
cither by using symmetrized (and normalized) wavefunctions to construct the density 
matrix or, equivalcntly, by averaging the distinguishable-particle density matrix over 
all N\ permutations (Feynman 1972): 
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Z= W\T, d N xp diat ({x u ...,x N },{x Pil) ,...,x PiN) },P) . (3.7) 
• P J 

This equation is illustrated in Fig. 13.41 for four particles. In each diagram, we arrange 
the points x\ , x% , X3 , X4 from left to right and indicate the permutation by lines. The 
final permutation (to the lower right of Fig. I3.4|) . corresponds for example to the 
permutation P(l,2, 3,4) = (4,3,2,1). It consists of two cycles of length 2, because 
P 2 (l) = P(4) = 1 and P 2 (2) = P(3) = 2. 
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Fig. 3.4 The 24 permutations contributing to the partition function for 4 ideal bosons. 

To illustrate eqn l|3.7[) . let us compute the contribution to the partition function 
stemming from this permutation for free particles: 

Z4321 = J dxidx 2 dx3dx4p llcc ({xi, x 2 , x 3 , X4} , {xa, x 3 , x 2 , xi} J3) = 
J dx 1 dx4p lroc (xi, X4, (3) p lroc (x4, x\, (3) 

f dx 1 p lla "(x 1 ,x 1 ,2f3)=z(2/3) sec cqn t3~2l 

This partial partition function of a four-particle system writes as the product of one- 
particle partition functions. The number of terms corresponds to the number of cycles 
and the length of cycles determines the effective inverse temperatures in the systems 
(Feynman 1972). 

The partition function for 4 bosons is given by a sum over 24 terms, 5 bosons involve 
120 terms, and the number of terms in the partition function for a million bosons has 
more than 5 million digits. Nevertheless, this sum over permutations can be computed 
exactly through an ingenious recursive procedure of N steps due to Landsberg (1961) 
(see Fig. UTS). 



J dx 2 dx 3 p hco (x 2 ,x 3 ,f3) p bcc {x 3 ,x 2 ,f3) 
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Fig. 3.5 The 24 permutations of Fig. 13.41 with the last-element cycle pulled out of 
the diagram. 
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Algorithm 3.2 canonic-recursion. py 



1 


import math 






2 


def z ( beta , k ) : 






3 


sum = 1/(1— math.exp(- 


-k* beta ) ) * *3 




4 


return sum 






5 


N=1000 






6 


beta = 1 ,/5 






7 


Z = [l.] 






8 


for M in range ( 1 ,N+ 1 ) : 






9 


Z . append (sum ( Z [ k ] * 


z(beta,Iv^k) for 


k in range (M))/M) 


10 


pi_list=[z(bcta ,k)*Z[N-k 


]/Z[N]/N for k 


in range ( 1 ,N+ 1 )] 



Figure 13.51 contains the same permutations as Fig. 13.41 but they have been rear- 
ranged and the last-clement cycle (the one containing the particle 4) has been pulled 
out. These pulled-out terms outside braces make up the partition function of a single 
particle at temperature (3 (on the first row), at inverse temperature 2/3 (on the second 
row), etc., as already computed in eqn (|3.6p . The diagrams within braces in Fig. 13.51 
contain the 3! = 6 three-boson partition functions making up Zj, (on the first line of 
the figure). The second row contains three times the diagrams making up Z 2 , etc. All 
these terms yield together the terms in eqn (|3.8I) : 

Zn = {z 1 Z N ^ 1 H h ZkZ N -k H h z N Z ) (ideal Bose gas) (3.8) 

(with Zq = 1, sec SMAC Sect. 4.2.3 for a proper derivation). Zq, as well as the single- 
particle partition functions are known from eqn (|3.6[) . so that we can first determine 
Z\, then i?2, and so on (see the Alg. harmonic-recursion. py, the SMAC web site 
contains a version with graphics output). 

The term in the Landsberg relation of eqn ()3.8|) oc ZfcZjv-fc can be interpreted as 
a cycle weight, the statistical weight of all permutations in which the particle N is in 
a cycle of length k 

^fc = TT^— z k z N-k (cycle weight). 

From the Landsberg recursion, we can explicitly compute cycle weights for arbitrary 
large ideal-boson systems at any temperature (sec Fig. I3.6|) . 

In Fig. 13.61 we notice that the cycle- weight distribution Ttk is flat for most k before 
it drops to zero around k ~ 780. Curiously, the derivative of this function yields the 
distribution of the condensate fraction (Holzmann and Krauth 1999, Chevallier and 
Krauth 2007), so that we see that at the temperature chosen, there are about 780 
particles in the groundstate (Nq/N ~ 0.78). At higher temperatures, the distribution 
of cycle weights is narrower. This means that there are no long cycles. 

We now turn the eqn (|3.8[) around, as we first did for the Gaussian integral: rather 
than computing Z(j3) from the cycle weights, we sample the cycle distribution from 
its weights that is, pick one of the 7Tfc with 

7Tl, . . . , 7Tfc, . . . , TT N . 

(we pick a cycle length 1 with probability 7Tx, 2 with probability 7T2, and so on (it is 
best to use the tower-sampling algorithm wc mentioned in the first lecture). Suppose 
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Fig. 3.6 Cycle-weight distribution for 1000 ideal bosons in a three-dimensional har- 
monic trap, obtained from Alg. canonic-recursion. py, exactly as written 

we sampled a cycle length k. We then know that our Bose gas contains a cycle of length 
k, and we can sample the three-dimensional positions of the k particles on this cycle 
from the three-dimensional version of Alg. harmonic-levy .py for k particles instead 
of N, and at inverse temperature k(3. Thereafter, we sample the next cycle length from 
the Landsberg recursion relation with N — k particles instead of N, and so on, until all 
particles are used up (see the Appendix for a complete program in 44 lines). Output 
of this program is shown in Fig. 13. 71 projected onto two dimensions. As in a real 
experiment, the three-dimensional harmonic potential confines the particles but, as 
we pass the Bose-Einstein transition temperature (roughly at the temperature of the 
left panel of Fig. 13. 7p . they start to move to the center and to produce the landmark 
peak in the spatial density. At the temperature of the right-side panel, roughly 80% 
of particles are condensed into the ground state. 



-6 6 




Fig. 3.7 Two-dimensional snapshots of 1000 ideal bosons in a three-dimensional 
harmonic trap (obtained with the direct-sampling algorithm). 

The power of the path-integral approach resides in the facility with which interac- 
tions can be included. This goes beyond what can be treated in an introductory lecture 
(see SMAC Sect. 3.4.2 for an in-depth discussion). For the time being we should take 
pride in our rudimentary sampling algorithm for ideal bosons, a true quantum Monte 
Carlo program in a nutshell. 
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Spin systems: samples and exact 
solutions 



4.1 Ising Markov-chains: local moves, cluster moves 

In this final lecture, we study models of discrete spins {<xi, . . . , ctjv} with = ±1 on 
a lattice with 2V sites, with energy 

E = ~J2 J ki°k<n- (4.1) 

(k,l) 

Each pair of neighboring sites k and I is counted only once. In eqn (|4.1[) . we may 
choose all the equal to +1. We then have the ferromagnetic Ising model. If we 
choose random values J^i = Jik — il 5 one speaks of the Edwards- Anderson spin 
glass model. Together with the hard-sphere systems, these models belong to the hall 
of fame of statistical mechanics, and have been the crystallization points for many 
developments in computational physics. Our goal in this lecture will be two-fold. We 
shall illustrate several algorithms for simulating Ising models and Ising spin glasses 
in a concrete setting. We shall also explore the relationship between the Monte Carlo 
sampling approach and analytic solutions, in this case the analytic solution for the 
two-dimensional Ising model initiated by Onsager (1942), Kac and Ward (1952), and 
Kaufmann (1949). 




Fig. 4.1 A spin flip in the Ising model. Configuration a has central field h = 2 (three 
"+" neighbors and one "— " neighbor), and configuration b has h = — 2. The statistical 
weights satisfy Trb/^a = exp (— 4/3). 



The simplest Monte Carlo algorithm for sampling the partition function 

Z= cxpH3£(<r)] 

confs (T 

picks a site at random, for example the central site on the square lattice of configuration 
a in Fig. 14.11 Flipping this spin would produce the configuration b. To satisfy detailed 
balance, eqn (|4.2|) . we must accept the move a —> b with probability min(l, irb/^a), as 
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Algorithm 4.1 markov-ising.py 



from random import uniform as ran , randint , choice 
from math import exp 

import square_neighbors # defines the neighbors of a site 
L = 32 
N = L*L 

S = [ choice ([ — 1 , 1] ) for k in range(N)] 
beta = 0.42 

nbr = square_neighbors (L) 
# nbr [k] = ( right , up , left , down ) 
for i.sweep in range (100): 
for iter in range(N): 
k = randint (0 ,N-1) 

h = sum (S [ nbr [k ] [ j ] ] for j in range(4)) 
Delta.E = 2.*h*S[k] 
Upsilon = exp(— beta * Delt a_E ) 
if ran(0.,l.) < Upsilon: S[k] =-S[k] 
print S 



implemented in the program Alg. markov-ising.py (see SMAC Sect. 5.2.1). This al- 
gorithm is very slow, especially near the phase transition between the low-temperature 
ferromagnetic phase and the high-temperature paramagnet. This slowdown is due to 
the fact that, close to the transitions, configurations with very different values of the 
total magnetization M = J2k ak contribute with considerable weight to the partition 
function (in two dimension, the distribution ranges practically from —N to N). One 
step of the local algorithm changes the total magnetization at most by a tiny amount, 
±2, and it thus takes approximately N 2 steps (as in a random walk) to go from one 
configuration to an independent one. This, in a nutshell, is the phenomenon of critical 
slowing down. 

One can overcome critical slowing down by flipping a whole cluster of spins simul- 
taneously, using moves that know about statistical mechanics (Wolff 1989). It is best 
to start from a random site, and then to repeatedly add to the cluster, with probability 
p, the neighbors of sites already present if the spins all have the same sign. At the end 
of the construction, the whole cluster is flipped. We now compute the value of p for 
which this algorithm has no rejections (see SMAC Sect. 5.2.3). 
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Fig. 4.2 A large cluster of like spins in the Ising model. The construction stops with 
the gray cluster in a with probability (1 — p) 14 , corresponding to the 14 "++" links 
across the boundary. The corresponding probability for b is (1 — p) 18 . 
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The cluster construction stops with the configuration a of Fig. l4~2l with probability 
p(a — > b) = const(l — p) 14 , one factor of (1 — p) for every link "++" across the 
cluster boundary. Likewise, the cluster construction in b stops as shown in b if the 

18 neighbors " " were considered without success (this happens with probability 

p(b — > a) = const(l — p) 18 ). The detailed balance condition relates the construction 
probabilities to the Boltzmann weights of the configurations 

TT a = const' exp [— /5(— 14 + 18)] 
n b = const' exp [—/?(— 18 + 14)] , 

where the const' describes all the contributions to the energy not coming from the 
boundary. We may enter stopping probabilities and Boltzmann weights into the de- 
tailed balance condition TT a p(a — > b) = 7Tf,7r(& — > a) and find 

exp (14/3) exp (-18/3) (1 - p) 14 = exp (-14/3) exp (18/3) (1 -p) 18 . (4.2) 

This is true for p = 1 — exp (—2/3), and is independent of our example, with its 14 

"++" and 18 neighbors " " links (sec SMAC Sect. 4.2.3). 

The cluster algorithm is implemented in a few lines (see Alg. cluster-ising.py) 
using the pocket approach of Section 12.31 Let the "pocket" comprise those sites of 
the cluster whose neighbors have not already been scrutinized. The algorithm starts 
by putting a random site both into the cluster and the pocket. One then takes a site 
out of the pocket, and adds neighboring sites (to the cluster and to the pocket) with 
probability p if their spins are the same, and if they are not already in the cluster. 
The construction ends when the pocket is empty. We then flip the cluster. 




Fig. 4.3 A large cluster with 1548 spins in a 64 x 64 Ising model with periodic 
boundary conditions. All the spins in the cluster flip together from "+" to "— " . 
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Algorithm 4.2 cluster-ising.py 



1 


from random import uniform as ran , 


randint , choice 


2 


from math import exp 




3 


import square neighbors # defines 


the neighbors of a site 


4 


L = 32 




5 


N = L*L 




6 


S = [ choice ([ — 1 , 1] ) for k in range 


(N)] 


7 


beta = 0.4407 




8 


p = 1— exp ( — 2* beta ) 




9 


nbr = square_neighbors (L) # nbr[k] 


= ( right , up , I eft , down ) 


10 


for iter in range(lOO): 




11 


k = randint (0 ,N-1) 




112 


Pocket = [k] 




13 


Cluster = [k] 




14 


while Pocket != [] : 




IS 


k = choice ( Pocket ) 




16 


for 1 in nbr [ k ] : 




17 


if S [ 1 ] = S [k] and 1 not 


in Cluster and ran (0,1) < p: 


18 


Pocket . append ( 1 ) 




19 


Cluster . append ( 1 ) 




20 


Pocket . remove ( k ) 




21 


for k in Cluster: S[k] =-S[k] 




22 


print S 





Cluster methods play a crucial role in computational statistical physics because, 
unlike local algorithms and unlike experiments, they do not suffer from critical slow- 
ing down. These methods have spread from the Ising model to many other fields of 
statistical physics. 

Today, the non-intuitive rules for the cluster construction are well understood, and 
the algorithms are very simple. In addition, the modern meta languages are so powerful 
that a rainy Les Houches afternoon provides ample time to implement the method, 
even for a complete non-expert in the field. 



4.2 Perfect sampling: semi-order and patches 

The local Metropolis algorithm picks a site at random and flips it with the probability 
min(l, lib/ Tr a ) (as in Section l4~Tj) . An alternative local Monte Carlo scheme is the heat- 
bath algorithm, where the spin is equilibrated in its local environment (see Fig. 14.41) . 
This means that in the presence of a molecular field h at site fc, the spin points up 
and down with probabilities itt and it7 , respectively, where 



2Sh ' 

(4.3) 



7T+ = = 

"h e -0E+ + e -0E- l_|_ e -2/3ft 



77 



e-^ E ~ 1_ 

e -f3E+ _|_ e -/3£- — i + g+2/3/i - 



The hcatbath algorithm (which is again much slower than the cluster algorithm, 
especially near the critical point) couples just like our trivial simulation in the five-site 
model of Section 11.21 At each step, we pick a random site k, and a random number 
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* + o 

* = * 

ifran(0,l)<7l^ 

o = o 

* = * 

else 
b 

Fig. 4.4 Heat bath algorithm for the Ising model. The new position of the spin k (in 
configuration b) is independent of its original position (in a). 

T = ran (0, 1), and apply the Monte Carlo update of Fig. 14. 41 with the same k and the 
same T to all the configurations of the Ising model or spin glass. After a time r coup , 
all input configurations yield identical output (Propp and Wilson 1996). 

For the Ising model (but not for the spin glass) it is very easy to compute the 
coupling time, because the half-order among spin configurations is preserved by the 
heat-bath dynamics (see Fig. 14. 5[) : we say that a configuration er = {(Ti, . . . , ctjv} 
is smaller than another configuration er' = {o~[, . . . ,cr' N } if for all k we have <7fc < 
crt. For the Ising model, the heat-bath algorithm preserves the half-order between 
configurations upon update because a configuration which is smaller than another one 
has a smaller field on all sites, thus a smaller value of ir£ . We just have to start the 
simulation from the all plus-polarized and the all minus-polarized configurations and 
wait until these two extremal configurations couple. This determines the coupling time 
for all 2 N configurations, exactly as in the earlier trivial example of Fig. 11.41 
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Fig. 4.5 Half-order in spin models: the configuration er_ is smaller and cr+ is larger 
than all other configurations, which are not all related among each other. 



The Ising spin glass does not allow such simple tricks to be played, and we must 
explicitly check the coupling for all 2 N initial configurations. This enormous surveying 
task can be simplified considerably by breaking up the configurations on the entire 
lattice into smaller "patches" (see Fig. I4.6[) . For a given patch size M (M = 16 in 
Fig. l4.6p . the total number of configurations on TV patches is bounded by N2 M <C 2 N . 



This is a half-order because not all configurations can be compared to each other. 
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We can now follow the heatbath simulation on all individual patches and at the end 
assemble configurations on the whole lattice in the same way as we assemble an entire 
puzzle picture from the individual pieces. The procedure can be made practical, and 
programmed very easily in modern meta languages such as Python. 

From a more fundamental point of view, the coupling concept in spin glass models 
is of interest because the physical understanding of these systems has seriously suffered 
from longstanding doubts about the quality of Monte Carlo calculations thus, at the 
most basic level, about the correct calculations of the correlation time. 





+ + - 


+ 




+ 


— + 


+ 


- + + 






+ + - 




+ 




+ 


+ 


— + 


+ 


+ 


- + + 






— + 





— + 

- + + 



— + - 



patch k 



spin configs 



patch / 



Fig. 4.6 The very large configuration space of an Ising spin glass on N sites viewed 
from the vantage point of N patches of smaller size M. 



4.3 Direct sampling, and the Onsager solution 

The two-dimensional Ising model is exactly solvable, as shown by Onsager (1942): 
we can compute its critical temperature, and many critical exponents. Let us be more 
precise: The partition function of the two-dimensional Ising model or the spin glasflon 
a planar lattice with N sites can be expressed as the square root of the determinant of 
a AN x AN matrix (see SMAC Sect. 5.1.4 for a practical algorithm). Periodic boundary 
conditions can also be handled, it gives four AN x AN matrices. This was used very 
successfully by Saul and Kardar (1992). For the Ising model on a finite square lattice, 
these matrices can be diagonalized analytically. This is the famous analytic expression 
for Z(/3) of Kaufman (1949) (see SMAC exerc. 5.9, p. 265). 

The solution of the Ising model amounts to summing its high-temperature series, 
and this solves an enumeration problem, as evidenced in the combinatorial solution 
by Kac and Ward (1949). Indeed, the density of states of the energy Af(E) can be 
extracted from the analytic solution (see Beale (1996)), this means that we can obtain 
the (integer) number of states for any energy E for two-dimensional Ising spin glasses 
of very large sizes. The exact solution of the Ising model thus performs an enumeration, 
but it counts configurations, it cannot list them. The final point we elaborate on in 
these lectures is that, unable to list configurations, we can still sample them. 

The subtle difference between counting and listing of sets implies that while we ac- 
cess without any trouble the density of states Af(E) we cannot obtain the distribution 
(the histogram) of magnetizations and know very little about the joint distribution of 

2 we speak here of one single "sample" of the spin glass, that is, a given choice of the couplings 
{Jkl}- To compute the average over all couplings is another matter. 
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energies and magnetizations, J\f(M,E). Doing so would in fact allow us to solve the 
Ising model in a magnetic field, which is not possible. 

Let us now expose a sampling algorithm (Chanal and Krauth 2009), which uses the 
analytic solution of the Ising model to compute M(M, E) (with only statistical, and 
no systematic errors) even though that is impossible to do (exactly). This algorithm 
constructs the sample one site after another. Let us suppose that the gray spins in the 
left panel of Fig. 14.71 are already fixed, as shown. We can now set a fictitious coupling 
J^i cither to — oo or two +oo and recalculate the partition function Z± with them. The 
statistical weight of all configurations in the original partition function with spin "+" 
is then given by 

= Z+exp(pj k i) 

+ Z+exp(PJki) + Z-exp(-pJkiy ' 
and this two- valued distribution can be sampled with one random number. Equa- 
tion (|4.4j) resembles the heatbath algorithm of eqn (|4.3p . but it is not part of a Markov 
chain: After obtaining the value of the spin on site k, we keep the fictitious coupling, 
and add more sites. Going over all sites, we can generate direct samples of the partition 
function of eqn (|4.ip . at any temperature, and with a fixed, temperature-independent 
effort, both for the two-dimensional Ising model and the Ising spin glass. This allows 
us to obtain arbitrary correlation functions, the generalized density of states Af(M, E), 
etc. We can also use this solution to determine the behavior of the Ising model in a 
magnetic field, which we cannot obtain from the analytic solution, although we use it 
in the algorithm. 
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Fig. 4.7 One iteration in the direct-sampling algorithm for the two-dimensional Ising 
model. The probabilities ir(<Tk = +1) and Tr(ak = — 1) (with k the central spin) are 
obtained from the exact solution of the Ising model with fictitious couplings J* k = 

±00. 



The correlation-free direct-sampling algorithm is primarily of theoretical interest, 
as it is restricted to two-dimensional Ising model and the Ising spin glass, which are 
already well understood. For example, it is known that the two-dimensional spin glass 
does not have a finite transition temperature. Nevertheless, it is intriguing that one 
can obtain, from the analytical solution, exactly and with a performance guarantee, 
quantities that the analytical solution cannot give. The "sampling" (even the very well 
controlled, full-performance-guarantee direct sampling used here) does not meet the 
limitations of complete enumerations. 
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Appendix A 



Here is an entire Quantum Monte Carlo program (in Python 2.5)) for ideal bosons in 
a harmonic trap (see the SMAC web site for a program with graphics output). 

Algorithm A.l direct-harmonic-bosons. py 



from math import sqrt , sinh , tanh , exp 
from random import uniform as ran , gauss 
def z ( beta , k ) : 

sum = 1/(1 — exp(—k* beta ))* *3 

return sum 
def cano ni c .recursion ( beta ,N ) : 

Z=[l.] 

for M in range ( 1 ,N+ 1 ) : 

Z . append (sum(Z [k] * z ( beta ,M-k) for k in range (M))/M) 
return Z 
def pi_list_make (Z,M) : 

pi.list =[0] + [z(bcta , k ) * Z pVHc ] / Z [M] /M for k in range (1 ,M+1)] 
pi_sum = [0] 

for k in range ( 1 ,M+ 1 ) : 

pi_sum . append (pi.sum [k — l]+pi_list [k]) 
return pi_sum 

def tower.sample ( data , Upsilon ) : #naive , cf . SMAC Sect. 1.2.3 
for k in range ( len ( data )) : 

if Upsilon<data [k ] : break 
return k 

def levy_harmonic_path ( Del_tau ,N) : # 
beta=N* Del_tau 

xN=gauss (0. , 1 . / sqrt (2* tanh ( beta /2. ) ) ) 
x=[xN] 

for k in rangc(l,N): 

Upsilon.l = l./tanh(Del_tau) + l./tanh((N-k)*Dcl_tau) 
Upsilon_2 = x [k-l]/sinh ( Dcl.t au )+xN/ s inh ((N-k)* Dcl.tau ) 
x_mean=Upsilon_2 /Upsilon_l 
sigma = l./ sqrt (Upsilon_l) 
x . append ( gauss ( x_mean , sigma ) ) 
return x 
N=512 

beta = l./2.4 

Z=c anon ic _r ec ur sion (beta ,N) 

x.config = [] 
y.config = [] 
while M > 0: 

pi_sum=p i _1 i st_make (Z ,M) 

Upsilon=ran (0,1) 

k=tower .sample ( pi_sum , Upsilon) 

x_config+=levy_harmonic_path (beta ,k) 

y _config+=levy_harmonic_path (beta ,k) 
M -= k 
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